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We study the dynamics of two-component atomic Bose gases initially in a mixture encountering a sudden 
quench of the inter-species interactions. The dynamics above the critical temperature Tc is studied using a lead- 
ing order large-N approximation that predicted a phase transition from mixing to phase separation as a function 
of the inter-species coupling. Here we explore the dynamics of this phase transition following a quench and 
compare our results to those found at zero temperature using the time-dependent Gross-Pitaevskii equations 
which ignore quantum and thermal fluctuations. In the regime above Tc where no condensate is present, how- 
ever, the time evolution of the densities following the quench exhibits features similar to that found at zero 
temperature where only the condensates contribute to the densities. When the inter-species interaction jumps 
above the critical value, we observe dynamical transitions from a homogeneous mixture to a phase-separated 
structure for both cases. Our simulations suggest that at temperatures above Tc where no condensate is present 
this dynamical transition should still be observable in experiments. 
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Experiments on tuning the interactions |[T] |2l and moni- 
toring the dynamics |3, 4) of two-component atomic Bose- 
Einstein condensates (BECs) have led to observations of in- 
teresting phenomena and inspired studies using such highly 
controllable systems. The dynamics of two-component Bose 
gases may simulate interesting phenomena in cosmology 
||5] |6l. In the weak-coupling regime at zero temperature 
(T = 0), two-component Bose gases have been broadly stud- 
ied using mean-field Bogoliubov theory showing a mixture to 
phase-separation transition when the inter-species interaction 
exceeds a critical value determined by the intra-species inter- 
actions ||7H3. 

The dynamics of multi-component BECs have been sim- 
ulated using the time-dependent Gross-Pitaevskii (TDGP) 
equations ifToHTSl . which are restricted to T and weak 
coupling due to the approximations involved. The stability 
prediction of the Bogoliubov theory (which is a one-loop cor- 
rection to the condensate contribution to the partition func- 
tion) is the same as that found by linear stability analysis of 
the TDGP equations. Self-consistent mean field theories allow 
one to study the dynamics at finite temperature as well as in 
stronger coupling regimes. Previously we have studied both 
the large-A^ expansion |[T3l and the leading-order-auxiliary- 
field (LOAF) theory fT4l for single-component Bose gases 
ifTSl . The LOAF theory becomes identical to the leading-order 
large-N expansion when there is no condensate. Thus the 
latter is appropriate for studying the dynamics of the phase- 
separation transition in the regime above the critical tempera- 
ture Tc, which has been less explored in the literature. 

A leading-order large-N approximation of two-component 
Bose gases shows that there is a mixture to phase-separation 
transition in equilibrium |T3l. Above Tc, there is no conden- 
sate and the transition is driven by a competition among the 
kinetic energy, intra-, and inter- species interactions. The goal 
of this paper is to study the dynamics of the system across this 
structural phase transition above Tc and compare it to the re- 
sults of a simulation at T = using the TDGP equations. We 
consider the case where an initial mixture of the two species 



evolves into a phase-separated structure after a sudden quench 
in the inter-species interactions. 

We first discuss zero temperature simulations of the quench 
using the TDGP equations to review what is know in that case, 
to benchmark our numerics, as well as to be able to compare 
with the finite temperature calculation. We will see below 
that the mode functions determining the coiTelation functions 
of the large-N approximation obey similar partial differential 
equations to those found for the condensates in the TDGP 
equations. We begin with the one-dimensional coupled TDGP 
equations [9| modeling a condensate at T = in a torus 

Here 0^, j = 1,2, denotes the condensate field for species 
j, i — 1 if i = 2 and j = 2 if j — 1, and the direction x 
is along the torus. We consider periodic boundary conditions 
and further assume Pi = P2 = Po ^nd mi = m2 = m. 
The unit of time is Iq = 2m/{hpQ^^) and the unit of length 

is ^^"^ = kg^. The coupling constants are related to the s- 
wave scattering lengths by A^j = ATrfi^aij/m for i,j — 1, 2. 
We will set h = 1. The TDGP equations can be solved using 
the split-step Fourier method |16, 17 1. We calibrate the grid 
size in real space and the time increment in our simulations 
using the exact solutions of Ref. flSl . 

For two-component BECs in equilibrium at T = 0, a struc- 
tural phase transition occurs when Af2/(AiiA22) > 1 |8, 9|. 
We consider the case where Af2/(AiiA22) < 1 initially so the 
system is a mixture. Then the inter-species interaction sud- 
denly changes to A^2/(AiiA22) > 1- Figure [l] shows the dy- 
namics of the density profiles and |02|^or (a) the case 
without any change of A12 and for (b) the case with a sudden 
quench in A12. We assume An and A22 are fixed and choose 
their values to be close to each other, which is an approxima- 
tion for the systems discussed in Refs. |1, 2|. In the initial 
density profiles we impose a small perturbation 1 19 1. Without 
a quench in A12 the profiles evolve smoothly with no observ- 
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Figure 1. (Color online) Evolution of the density profiles for two- 
component BEC according to the TDGP equations ([T]l. The ini- 
tial state is at r = and fcoan = 0.11, fcoa22 = 0.1. The 
profiles are offset by 2 for each time increment and we show t = 
0,26.2fo,52.4to. (a) A12/A22 = 0.5. (b) A12/A22 = 1.5. The 
structure of (b) is shown in more detail in (c) for t — 52.4to without 
offset. The solid (dashed) line correspond to species 1 (2). 



able structure. When A12 suddenly jumps above the critical 
value, the two components start to form domains and exclude 
each other as shown in Fig. [T](b) and (c). The large number 
of domains emerging during the dynamics suggests that the 
initial mixture state is a highly excited state when compared 
to the genuine ground state of the large value of Ai2. Ref. ifTIl 
investigated more dynamics using TDGP equations. 

In the BEC phase, to go beyond the zero-temperature and 
weak-coupling limit of the TDGP equations, one should in- 
clude the pair correlations discussed in Refs. fT4"20l as well 
as the fluctuations of the normal density included in the large- 
N approximation. The interplay among the condensate, quasi- 
particles, and pair correlations during the dynamics will be an 
interesting topic for future studies. 

In the finite-temperature quantum regime above Tc, all con- 
densates vanish and only the fluctuations of the normal density 
are important. The relevant mean-field theory which preserves 
symmetries is the leading term in the large-N expansion of the 
theory, which for the normal phase of a two-component Bose 



gas was presented in Ref. (13). There we focused on the static 
properties of the theory and the partition function which was 
relevant to study the phase diagram. Here we study the dy- 
namics coming from the effective action derived in Ref. [ 13 1. 
The action describing two-component bosons with contact in- 
teractions is 



5 = y" [dx] 
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Here [dx] = dtd'^x. We introduce the large parame- 
ter N into the theory by the replication trick (f>j,n 
where n = 1,2, ...N and rescale the coupling constants 
Xjj — >■ Xjj/N and A12 — >■ ^12/^ lfT3l . The generating 
functional for the correlation functions is given by Z[J] = 

I {Uf=i,2-n=i^'l^j,n'D'l^ln) e''^!^'"^-"'^.* "1 . Here the action 
with the source term after the replication becomes 
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Here <i> = {(jii 

dtagih[+\h[-\hi+\hi-'> 
diag{fj,i, ^i, ^2, f^2, ■ ■ ■) is the bare Green's function, 
and hf ^ = Tidt - V'^/{2mj) for j = 1,2. There are N 

copies in G(7^, and Gq ^. J is the source coupled to 
We next insert the identity 



l[VxPa, I e 
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inside the path integral for Z[J] which introduces two func- 
tional delta functions fTPi so one can replace X)!r=i 't'j n4>j-n 
by {N/Xjj)aj in S. Here the x integration contour runs par- 
allel to the imaginary axis (see Ref. lUTj). This replacement 
facilitates our resummation scheme and we will treat as a 
small parameter Let G^^ = Gq^ +diag{xi, Xij X2, X2, ■■ ■)■ 
After performing the Gaussian integral in the (j>j^n^ one has 

Seff = 
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Here X = (xi,X2,ai,a2)"^ with its source term K and ex- 
pectation value Xc- We evaluate the path integrals over Xj, ctj 
via the method of stationary phase and at leading order we 
keep only the contributions at the stationary phase point. 
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The generating functional of the one-particle irreducible di- 
agrams is obtained from In Z by Legendre transform: F — 
l\dx\ ( J'f <i>c + K'^Xc + i In Z), where $c and are the clas- 
sical values of $ and X respectively. Keeping the leading 
term in the expansion and then setting N = 1, we obtain 



r = / {dx\ 
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Here $ = ( 
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and Gn has been reduced to a 



4x4 matrix. The equation for the condensate in the presence 
of sources is given by (5r/(5$^^ = J = J[da:]G(^^$. 

For the initial conditions of our simulations we consider 
static homogeneous fields and define the effective potential 
as Veff — T/{NVt), where i is the length of time in the 
integral. The broken-symmetry condition is determined from 
the true minimum of the effective potential: 5Vef f /Scj)* — 0, 



which becomes Xj4>j = 0- In the normal phase 



while 



in the broken-symmetry phase Xj — 0. In the normal phase, 
the term ^^^Gq ^<1> is zero at the minimum of the potential 
(which occurs at (pj — 0). The condition 5Veff/Saj = 
fixes the relation between aj and fij . Thus 
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Here the expectation value means Tr[pi)(f>*(j)j ], where po is 
the density matrix describing the quantum state at the initial 
time. The parameters are determined self-consistently at the 
initial time (see below). In the BEC phase the expectation 
value would include the condensates plus the quantum fluctu- 
ations. 

To describe the dynamics, we derive the equations of mo- 
tion from the identity j dzG^^{x,z)G{z,y) = (5c(x — y), 
where G^^ = S'^T /S(paS4'b and we use the closed time path 
(CTP) boundary conditions on the Green's functions due to 
Schwinger and Keldysh (S-K) (see ll22l |23 | and references 
therein). Here (j)a,b = 'Pjy ^c{^ — y) is a four-dimensional 
delta function defined on the S-K contour, and x, y,z de- 



note four-vectors. One finds that Gj^^{\,y) 



y){h^^^ + Xj) and the off-diagonal elements vanish. Fol- 
lowing Ref. f22], the equations for the Green's functions can 
be satisfied by introducing a set of fields (f>j{x) satisfying 
the homogeneous equation J dyGj^^ {\,y)(j)j{y) — 0, where 
dy = dtd'^y. Explicitly, 

idt(t>j = {-dl - Hj+c 
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Here aj = Xjj 



'j(j>j) and the expectation value is over the 
initial (thermal) state of the system. This is to be compared 
with Eq. ([T]) for the T = case. 

To implement the initial condition on the Green's function, 
one represents the quantum field (f>i and its conjugate by a 
mode expansion in terms of canonical annihilation and cre- 
ation operators: 0j(x) = ^^-^ a t). The operator 
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Figure 2. (Color online) Evolution of the density profiles for species 
1 (solid lines) and 2 (dashed lines). Here koau = 0.1, koa22 = 
0.11, r = 3To, and iV = 256. Initially koai2{t = 0) = 0.5. The 
profiles for t > is offset by 2 after each time increment. We show 
the profiles at t = 0, 5.2to, and lO.bto. (a) shows A;oai2(t > 0) = 
0.5 below the critical value, (b) shows fcoai2(i > 0) = 1.5 above 
the critical value. The detailed structure of (b) is shown in (c) and 
the black line shows the density difference Ap — pi — p2. 



Ujk is the time independent initial value of the annihilation 
operator in the Heisenberg picture. The temperature enters 
by specifying that the initial distribution for each species is 
given by a Bose-Einstein distiibution at a given temperature. 

,t \ _ x...x...i^..^f..../i^^^\ ii-i 



Thus {a'^k^-j'k') = %'4fc'[cxp(wjo/fcBr) - l]-\ where 
ujjo = k'^/{2mj) + Xj{t = 0). In our simulations we fo- 
cus on the motion along the direction along the torus and ne- 
glect the details in the other two spatial dimensions. The time- 
ordered Green's function of the j-th species is then written in 
terms of these fields as Gjj(x,y) = —i{Tc[4>j{'^)(t>j{y)])^ or 
explicitly ^^(xo - yo)G'>j(x,y) + Ociyo - xo)G<j(x,y). 
Here 7c is the CTP time-ordered product and 9c {^0) is the 
CTP step function while G>.j(x,y) — —i{ipj{x)(j>*{y)) and 
G<j = (x)(/)j (y)) following the convention of Ref. |23|. 
In terms of the fields (f>j{x.), we have Pj{x) = (0*(x)0j(x)). 

In our simulations, we follow N modes, where TV is also 
the number of grid points in x. The initial conditions are 
chose as fjk{x,0) = y^T/Nj^ exp[i{kx — i^jot)], where the 
normalization coefficients are chosen to satisfy Pj{x,t = 
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0) = \^jix,t = 0)|2 = po = fco for j = 1,2. The 
equations of motion (|7]i then lead to idtfjk{x, t) — [—d1 + 
Xj{x,t)]fjk{x,t). The densities are evaluated by (a;, i) — 
{^*{x,t)dpj{x,t)) = Y.kfU^^i)fok{x,t)n,k{T). The com- 
posite fields are then evaluated by aj{x, t) = XjjPj{x, t) and 
Xj{x,t) = —^j+aj{x,t) + ^a-j{x,t). Here we assume the 

relation between aj (x, t) and Xj {x, t) follows the equilibrium 
relation (|6]l but the values of /ij are fixed by the equilibrium 
values so that the evolution is continuous. 

We again consider the setup where initially the system is 
a mixture described by a thermal distribution at temperature 
T > Tc and suddenly A12 jumps above the critical value. 
The initial condition is determined by solving the equilibrium 
equations of state pj = —dVef f /dpj and Eq. (|6]l tlSl: 

Pj = = 0) = ^jjPj- (8) 

p 

Here n(a;) = [exp(a:;/A:BT) — ,ujj = p"^ /{2mj)+Xj(t = 
0), and p is a continuous variable representing the momen- 
tum. We choose pi = /92 = Po ™d nii = 1122 ~ m. 
Eq. (|8]l gives the values for Xj ™d aj, and /ij can be de- 
termined from Eq. (|6]|. For fcooii = 0.11, A;oi22 — 0.1, and 
^0^12 (i = 0) = 0.5, the system is stable in the mixture phase 
and the corresponding parameters are obtained. For illustra- 
tive purposes, the initial temperature is chosen as T = 3To, 
where Tq is the BEC temperature of an noninteracting Bose 
gas with density pc |24|. The dynamics when fcoai2(< > 0) 
jumps to 1.5 (above the critical value) is then simulated. 

The evolution of the density profiles is shown in Figures |2] 
(a) and (b) corresponding to the cases without and with a jump 
in A12 for N = 256 modes. In contrast to the T = Q case, here 
we can use a uniform initial condition with pi = P2- One can 
see that if A12 remains below the critical value (which is close 



to 10\/AiiA22), the system remains a mixture. In contrast, 
when A12 jumps above the critical value, the two density pro- 
files evolve into domains and the domains of different species 
avoid each other Therefore many interfaces separating differ- 
ent species emerge during the dynamics. We found that our 
results are robust when the number of modes double or when 
smaller time steps are taken. For larger jumps of A12 we found 
that the time it takes to form those domains decreases. If T/Tq 
is too large, the kinetic energy dominates and one cannot find 
a structural transition for reasonable parameters. When com- 
pared to the results at T = 0, one can see that the modula- 
tions of the densities at finite temperature above Tq are not as 
strong. Therefore above To we expect that density ripples cor- 
responding to the phase-separated phase should emerge when 
A12 is above the critical value. The mixture background above 
To is due to kinetic energy effects that are not significant at 
T = 0. The density difference shown on Fig.|2]should be large 
enough to be distinguished from the mixture background in 
experiments. For larger jumps of A12, the density ripples be- 
comes more prominent due to larger inter-species interactions 
so experimentally one should use the largest available jump 
for better observations of this dynamical structural transition. 

In summary, we found that a dynamical structural transition 
should be observable in two-component Bose gases at T = 
and temperatures above Tc so the presence of BECs is not a 
crucial ingredient for this transition. Our formalism for finite 
T dynamics provides a coherent description and complements 
conventional kinetic -theory approaches lIZSll . If optically con- 
trollable collisions (see references in Ref. ^2611) can be realized 
in two-component Bose gases, one may study more compli- 
cated structures and dynamics in the presence of inhomoge- 
neous interactions. 
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